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Abstract 

Hysteresis  and  constitutive  nonlinearities  are  inherent  properties  of  ferroelectric  transducer  ma¬ 
terials  due  to  the  noncentrosymmetric  nature  of  the  compounds.  In  certain  regimes,  these  effects 
can  be  mitigated  through  restricted  input  fields,  charge-  or  current-controlled  amplifiers,  or  feed¬ 
back  designs.  For  general  operating  conditions,  however,  these  properties  must  be  accommodated 
in  models,  transducer  designs,  and  model-based  control  algorithms  to  achieve  the  novel  capabilities 
provided  by  the  compounds.  In  this  paper,  we  illustrate  the  construction  of  inverse  filters,  based  on 
homogenized  energy  models,  which  can  be  used  to  approximately  linearize  the  piezoceranric  trans¬ 
ducer  behavior  for  linear  design  and  control  implementation.  Attributes  of  the  inverse  filters  are 
illustrated  through  numerical  examples  and  experimental  open  loop  control  implementation  for  an 
atomic  force  microscope  stage. 
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1  Introduction 


Ferroelectric  materials,  including  the  compound  lead  zirconate  titanate  (PZT),  exhibit  novel  actuator 
and  sensor  capabilities  due  to  the  unique  electromechanical  coupling  which  they  exhibit.  This  pro¬ 
vides  them  with  the  capability  for  providing  broadband  transduction  and  nanometer-level  set  point 
accuracy.  Furthermore,  piezoelectric  transducers  are  moderately  inexpensive  and  can  be  designed  to 
minimally  affect  the  passive  dynamics  of  underlying  structures.  However,  the  noncentrosymmetric 
ion  structure  that  imbues  the  materials  with  unique  actuator  and  sensor  properties  also  produces 
hysteresis  and  constitutive  nonlinearities  at  all  drive  levels. 

To  illustrate,  consider  the  prototypical  atomic  force  microscope  (AFM)  stage  depicted  in  Fig¬ 
ure  1(a)  which  employs  stacked  piezoceramic  actuators  to  position  the  sample  in  the  x  and  y  di¬ 
rections.  An  additional  PZT  mechanism  provides  transverse  positioning  capabilities.  Nested  minor 
loops  collected  at  0.1  Hz  are  plotted  in  Figure  1(b)  and  data  collected  at  frequencies  ranging  from 
0.279  Hz  to  27.9  Hz  is  plotted  in  Figure  2  to  illustrate  the  frequency-dependent  nature  of  the  hys¬ 
teresis  inherent  to  field-displacement  data. 

At  low  frequencies,  the  inherent  hysteresis  can  be  accommodated  through  proportional-integral- 
derivative  (PID)  or  robust  control  designs  [5,6,15,19].  However,  at  the  higher  frequencies  required 
for  applications  ranging  from  real-time  monitoring  of  biological  processes  —  e.g.,  protein  unfolding 
—  to  comprehensive  product  diagnostics,  increasing  noise-to-data  ratios  and  diminishing  high-pass 
characteristics  of  control  filters  preclude  a  sole  reliance  on  feedback  laws  to  eliminate  hysteresis. 

Alternatively,  it  is  illustrated  in  [12,13],  that  use  of  charge-  or  current-controlled  amplifiers  can 
essentially  eliminate  hysteresis.  However,  this  mode  of  operation  can  be  prohibitively  expensive  when 
compared  with  the  more  commonly  employed  voltage-controlled  amplifiers,  and  current  control  is 
ineffective  if  maintaining  DC  offsets  as  is  the  case  when  the  x-stage  of  an  AFM  is  held  in  a  fixed 
position  while  a  sweep  is  performed  with  the  y-stage. 

This  motivates  the  development  of  models  and  model-based  control  designs  which  incorporate 
and  accommodate  the  hysteresis  and  constitutive  nonlinearities.  Numerous  approaches  have  been 
employed  to  characterize  these  nonlinear  effects  including  Preisach  models  [7, 18] ,  domain  wall  mod¬ 
els  [25,26],  micromechanical  models  [4,10,11],  mesoscopic  energy  relations  [3,9]  and  homogenized 
energy  models  [23,30].  We  employ  the  homogenized  energy  framework  due  to  its  energy  basis,  its 
capability  to  quantify  a  wide  range  of  physical  phenomena  and  operating  regimes,  its  unified  nature 


Figure  1:  (a)  PZT-based  AFM  stage,  and  (b)  nested  minor  loops  in  data  collected  at  0.1  Hz. 
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Figure  2:  Frequency-dependent  field-displacement  behavior  of  a  stacked  PZT  stage:  sample  rates  of 
(a)  0.279  Hz,  (b)  5.58  Hz,  and  (c)  27.9  Hz. 


for  characterizing  hysteresis  in  ferroelectric,  ferromagnetic  and  ferroelastic  compounds  [28,29],  and 
the  potential  it  provides  for  real-time  implementation.  Details  regarding  the  development  of  this 
modeling  framework  and  its  relation  to  other  characterization  techniques  can  be  found  in  [21,30]. 

Model-based  control  design  for  piezoceramic  transducers  operating  in  highly  hysteretic  and  non¬ 
linear  regimes  can  be  roughly  segregated  into  two  categories:  (i)  nonlinear  control  designs,  and  (ii) 
linear  control  designs  employing  nonlinear  inverse  filters.  Examples  of  the  first  technique  in  the 
context  of  optimal  control  design  for  smart  material  transducers  are  provided  in  [20,32].  The  second 
technique  is  based  on  the  concept  of  employing  either  an  exact  or  approximate  inverse  model  to  lin¬ 
earize  the  transducer  behavior  in  the  manner  depicted  in  Figure  3.  This  approach  has  been  employed 
with  a  variety  of  models  and  control  designs  —  e.g.,  see  [31]  for  details  regarding  the  development 
of  adaptive  control  designs  utilizing  piecewise  linear  Preisach  models  and  their  inverses  —  and  is  the 
technique  which  we  focus  on  in  this  paper. 

In  Section  2  we  summarize  constitutive  relations  developed  in  [21,  30]  for  ferroelectric  materi¬ 
als  and  provide  a  highly  efficient  algorithm  for  implementing  the  model  when  thermal  relaxation 
is  negligible.  A  corresponding  inverse  polarization-field  algorithm  is  summarized  in  Section  3  and 
illustrated  through  a  numerical  example.  The  constitutive  model  is  subsequently  employed  in  Sec¬ 
tion  4  to  develop  a  lumped  model  for  the  stacked  actuator  employed  in  the  AFM  stage  shown  in 
Figure  1(a)  to  illustrate  the  construction  of  a  macroscopic  transducer  model.  The  accuracy  of  the 
transducer  model  is  illustrated  through  a  comparison  with  the  frequency-dependent  data  plotted 
in  Figure  2.  In  Section  5,  an  algorithm  for  the  inverse  displacement-field  relation  to  linearize  the 
transducer  response  is  developed  and,  in  Section  6,  the  algorithm  is  experimentally  implemented 
as  an  inverse  filter  for  the  open  loop  tracking  of  a  triangular  input  signal.  It  is  demonstrated  that 
this  model-based  filter  design  effectively  linearizes  the  nonlinear  and  hysteretic  transducer  dynamics 
and  provides  an  approximately  tenfold  increase  in  accuracy  at  higher  frequencies  as  compared  with 
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Figure  3:  Use  of  an  inverse  filter  to  linearize  the  response  a  of  a  hysteretic  actuator  to  achieve  a 
desired  output  Ud- 
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the  unfiltered  case.  This  significantly  improves  the  accuracy  of  the  transducer  and  diminishes  the 
sole  reliance  on  feedback  laws  whose  authority  decrease  as  tracking  speeds  and  noise-to-data  ratios 
increase. 


2  Constitutive  Relations 

Constitutive  relations  quantifying  the  electromechanical  behavior  of  piezoceranric  materials  are  con¬ 
structed  in  two  steps.  In  the  first,  Helmholtz  and  Gibbs  energy  relations  at  the  lattice  level  are  used 
to  characterize  the  local  field-polarization  and  field-strain  behavior  of  ferroelectric  compounds  for 
thermally  inactive  and  active  operating  regimes.  In  the  second  step  of  the  development,  material 
nonhomogeneities  and  variable  effective  field  effects  are  incorporated  through  the  assumption  that 
certain  material  properties  are  manifestations  of  underlying  distributions  rather  than  constants.  This 
yields  low-order  macroscopic  constitutive  relations  which  are  efficient  to  implement. 

2.1  Local  Constitutive  Relations 

Let  E ,  P,  e  and  a  respectively  denote  the  electric  field,  polarization,  strain  and  stress.  It  is  illustrated 
in  [30]  that  an  appropriate  formulation  of  the  Helmholtz  energy  for  fixed  temperatures  in  the  absence 
of  stresses  is 

'  \n{P  +  PR )2  ,  P< -Pi 

il>P(P)  =  <  h(p  ~  pR )2  ,  P>  Pi 

k  Hpi  -  pR)  (S  -  Pr)  >  1^1  <  Pi- 

As  shown  in  Figure  4,  Pj  is  the  positive  inflection  point  which  delineates  the  transition  between 
stable  and  unstable  regions,  Pq  denotes  the  unstable  equilibrium,  and  Pr  is  the  value  of  P  at  which 
the  positive  local  minimum  of  i/)  occurs.  The  parameter  r/  is  the  reciprocal  of  the  slope  of  the  E-P 
relation  after  switching  occurs.  This  fact  can  be  used  to  establish  an  initial  parameter  value  for  r/ 
when  modeling  a  specific  data  set. 


V(P)=G(0,P) 


Figure  4:  (a)  Helmholtz  energy  if)  and  Gibbs  energy  G  for  a  =  0  and  increasing  fields  E.  (b)  Switch 
in  the  local  polarization  P  that  occurs  as  E  is  increased  beyond  the  local  coercive  field  Ec  given  by 
(5)  in  the  absence  of  thermal  activation. 
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The  corresponding  Gibbs  energy  relation 

GP(E,P)  =  ^p(P)-EP  (1) 

incorporates  the  electrostatic  energy  due  to  the  applied  field  E  when  a  =  0. 

Elastic  effects  and  electromechanical  coupling  are  incorporated  in  the  Helmholtz  energy  relation 

ip(P,  e)  =  ipp(P)  +  ^Ye2  -  clieP  -  a2eP2 ■ 

The  Gibbs  energy  is  then  given  by 

G(E,  a,  P ,  e)  =  ipP{P)  +  ^ Ye 2  -  oi eP  -  a2eP 2  -  EP  -  ae  (2) 

where  as  incorporates  the  elastic  energy.  Note  that  Y  denotes  the  Young’s  modulus  and  a\,a2  are 
ferroelastic  coupling  coefficients  associated  with  linear  piezoelectric  and  quadratic  electrostrictive 
effects. 


Polarization  Kernel  —  Negligible  Thermal  Activation 

In  the  case  of  negligible  thermal  activation,  the  local  average  polarization  kernel  P  is  determined 
from  the  necessary  conditions 

dG  n  d2G  n 
dP  ~  ’  5P2  >  • 

Applying  these  conditions  to  (1)  yields  the  piecewise  linear  E-P  relation 


P(E)  =  -E  +  Pr5  (3) 

V 

where  5  =  —  1  for  negatively  oriented  dipoles  and  (5  =  1  for  those  with  positive  orientation.  To  specify 
<5,  and  hence  P,  more  specifically  in  terms  of  the  initial  dipole  orientations  and  previous  switches, 
we  employ  Preisach  notation  and  take 


[p(E-,Ec,om 


'  [P(E-,Ec,0m  ,  r(t)  =  0 

<  f  ~  Pr  ,  r(t)  /  0  and  £(maxr(t))  =  -Ec 

.  f  +  Pr  ,  r(t)  /  0  and  £(maxr(t))  =  Ec. 


(4) 


Here 


[P(E;Ec,Om 


l  f  +  Pr 

defines  initial  kernel  values  in  terms  of  the  parameter  £ 
the  set  of  switching  times  is  given  by 


,  E(0)  <  —Ec 
,  -Ec  <  E( 0)  <  Ec 
,  E( 0)  >  Ec 

=  rj-  ±  Pr,  0  designates  the  empty  set,  and 


r(t)  =  {ts  6  (0,  i]  |  E(ts)  =  -Ec  or  E(ts)  =  Ec}. 


The  local  coercive  field 

Ec  =  v(Pr  ~  Pi)  (5) 

quantifies  the  field  at  which  the  negative  well  ceases  to  exist  and  hence  a  dipole  switch  occurs.  To 
illustrate,  the  condition  r  ^  0  and  E(maxr(f))  =  Ec  designates  that  switching  has  occurred  and  the 
last  switch  was  at  Ec ;  hence  the  local  polarization  is  \P(E\  Ec,  £)](£)  =  +  Pr 
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Polarization  Kernel  —  Thermal  Activation 


If  thermal  activation  is  significant,  dipoles  can  achieve  the  thermal  energy  required  to  switch  in 
advance  of  the  minimum  Gibbs  energy  so  the  relative  thermal  energy  kT /V  and  Gibbs  energy  G 
must  be  balanced  through  Boltzmann  principles.  The  probability  density  for  achieving  an  energy 
level  G  is  given  by 

p{G)  =  Ce~GV/kT  (6) 

where  k  is  Boltzmann’s  constant,  V  is  a  reference  volume  and  C  is  a  constant  that  is  selected  so  that 
when  n(G)  is  integrated  over  all  possible  dipole  orientations,  a  probability  of  unity  is  achieved.  If  we 
let  2e  be  the  separation  between  possible  polarization  states  around  Po ,  the  probabilities  of  reaching 
a  polarization  state  having  sufficient  energy  to  switch  orientations  are  given  by 

rpO+e  -G(E,P)V/kTdp  fPo+e  -G(E,P)V/kT  dp 

T+~  ~  /~_£  e~G(E,P)V/kTdp  ’  r-+  “  jPo+e  e-G(E,P)V/kTdp  '  l  j 


The  likelihoods  of  reaching  the  required  energy  and  thus  of  the  dipoles  switching  from  a  positive  to 
a  negative  orientation  and  conversely  are  then 


(8) 


where  T(T)  is  the  relaxation  time  at  temperature  T.  The  fractions  of  dipoles  in  each  orientation 
evolve  according  to  the  ordinary  differential  equations 

dx+ 

—  =  -p+-x+  +p-+x- 
dx _ 

—  =  -P-+X-  +P+-X+. 

The  expected  polarizations  due  to  positively  and  negatively  oriented  dipoles  are 


POO 

(P+)  =  /  Pp(G)dP 

■'  Po+e 

so  that  evaluation  of  C  yields 


(P+) 


Jpo+ePe~G^p^v/kTdP 

fh  e~G{E^T)VlkTdp 

JPo+e 


The  local  average  polarization  is  subsequently 


(p+) 


Pp(G)dP 


(P-) 


fpo-c  pe~G(E,P,T)V/kT dp 
J  — oo 

rp [)-£  e-G(E,P,T)V/kTdp 
J— OO 


(9) 


P  =  x+{P+)  +  x-(P-).  (10) 

In  the  manner  detailed  in  [30],  the  evaluation  of  the  integrals  in  (7)  and  (9)  can  be  simplified  through 
approximations  employing  the  inflection  points  PPj  rather  than  the  unstable  equilibrium  Po- 

2.2  Global  Constitutive  Relations 

For  homogeneous  compounds  with  uniform  effective  fields  Ee,  the  local  lattice  relations  (3),  (4) 
or  (10)  can  be  extrapolated  throughout  the  material  to  provide  global  constitutive  relations.  This 
yields  the  nearly  instantaneous  transitions  at  coercivity  that  are  associated  with  certain  single  crystal 
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compounds  —  e.g.,  the  hysteresis  kernels  depicted  in  Figure  4  provide  a  reasonable  characterization  of 
the  single  crystal  BaTiC>3  behavior  shown  on  pages  72-76  of  [14]  —  but  provide  a  poor  characterization 
of  the  mollified  transition  behavior  of  general  ferroelectric  compounds.  To  incorporate  the  effects  of 
material  nonhomogeneities,  poly  crystallinity,  and  variable  effective  fields  Ee  =  E  +  Ej,  we  assume 
that  the  interaction  field  Ej.  due  to  neighboring  dipoles  and  certain  electromechanical  interactions  [1], 
and  local  coercive  field  Ec  given  by  (5)  are  manifestations  of  underlying  distributions  rather  than 
constants.  If  we  designate  the  associated  densities  by  v\  and  v 2,  the  macroscopic  field-polarization 
behavior  is  quantified  by  the  relation 

roo  /*oo 

[P(E)](t)=  /  /  [P(E  +  (11) 

J  0  J  —  00 


where  the  kernel  P  is  given  by  (3),  (4)  or  (10). 

As  detailed  in  [23],  the  densities  v\  and  v2  are  assumed  to  satisfy  the  physical  criteria 

(i)  v\ (x)  defined  for  x  >  0, 

(ii)  v2(-x)  =  v  2(2), 

(hi)  Mx)|  <  cie~aix, 

h(x)j  <  c2e~a2^ 


(12) 


for  positive  ci,  01,02,02-  The  restricted  domain  in  (i)  reflects  the  fact  that  the  coercive  field  Ec  is 
positive  whereas  the  symmetry  enforced  in  the  interaction  field  through  (ii)  yields  the  symmetry 
observed  in  low- field  Rayleigh  loops.  Hypothesis  (iii)  incorporates  the  physical  observation  that  the 
coercive  and  interaction  fields  decay  as  a  function  of  distance  and  guarantees  that  integration  against 
the  piecewise  linear  kernel  yields  finite  polarization  values. 

By  employing  numerical  integration  routines  tailored  to  the  infinite  domains  or  truncated  intervals 
resulting  from  the  decay  criteria  (12),  the  integrals  in  (11)  can  be  approximated  to  obtain  the 
discretized  model 

Ni  Nj 

[P(E)](t)  =  ^^\P(E  +  EI.-,ECi^j)]^(E<H)u2(EI.)viwj.  (13) 

*= 1  3= 1 


Specific  choices  for  the  weights  Vi,Wj  and  abscissas  Eej,ECi  are  detailed  in  [21,30]. 

Techniques  for  identifying  the  densities  v\  and  z^2  are  illustrated  in  [23].  For  certain  applications, 
reasonable  accuracy  is  provided  by  a  priori  functions  satisfying  the  physical  criteria  (12)  and  having 
a  small  number  of  parameters  to  be  estimated  through  least  squares  fits  to  data  —  e.g.,  variances 
and  means  in  normal  and  lognormal  relations.  For  more  general  applications  requiring  high  accuracy 
for  a  wide  range  of  operating  conditions,  the  Ni  +  Nj  discretized  density  values  v\  (ECi)  and  u2 (£/-) 
can  be  estimated  through  least  squares  techniques. 

To  obtain  an  elastic  constitutive  relation,  the  equilibrium  condition 


is  invoked  to  obtain 

a  =  Ye  —  a\P  —  a2P2.  (14) 

When  P  =  0,  (14)  reduces  to  Hooke’s  law.  To  incorporate  internal  damping,  we  posit  that  when 
P  =  0,  stress  is  proportional  to  a  linear  combination  of  strain  and  strain  rate  (Kelvin-Voigt  damping 
hypothesis).  This  yields  the  constitutive  relation 

a  =  Ys  —  cd£  —  a\P  —  a2P 2  (15) 
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where  cd  is  the  Kelvin-Voigt  damping  parameter.  The  combination  of  the  field-polarization  model 
(11)  or  (13)  and  the  electromechanical  relation  (15)  are  employed  in  Section  4  to  construct  a  lumped 
model  for  a  stacked  PZT  actuator  operating  in  hysteretic  and  nonlinear  regimes. 


2.3  Implementation  Algorithm 

The  efficiency  of  inverse  algorithms  used  to  construct  inverse  filters  is  dependent  on  the  efficiency  of 
forward  algorithms  used  to  implement  the  discretized  polarization  model  (13).  In  this  section,  we 
summarize  a  highly  efficient  algorithm  to  evaluate  (13)  when  the  kernel  P  is  given  by  (3)  or  (4)  for 
regimes  in  which  relaxation  processes  are  negligible.  Analogous  algorithms  for  the  thermally  active 
kernel  (10)  are  reported  in  [2]. 

From  (4),  it  is  observed  that  for  each  field  value  Ej:j ,  it  is  necessary  to  determine  whether  a 
transition  has  occurred  relative  to  the  coercive  value  ECi.  This  yields  Ni  x  Nj  conditions  to  be 
checked  for  each  input  value.  While  this  can  be  easily  accomplished  using  an  if-then  construct, 
implementation  in  this  manner  diminishes  significantly  the  efficiency  of  the  algorithm.  This  motivates 
consideration  of  an  algebraic  technique  for  evaluating  the  conditional  statements. 

To  retain  the  history  of  whether  or  not  effective  field  values  Eej  =  E  +  Ejj  have  switched  due  to 
encounters  with  coercive  field  values  Ec% ,  we  employ  (3)  to  motivate  the  matrix  formulation 

P=-  +  PrA(E;Ec,Ei) 

T] 

where  A  =  1  if  evaluating  on  the  upper  branch  of  the  hysteron  (hystersis  kernel)  and  A  =  —1  if  on 
the  lower  branch.  For  the  evaluation  of  (13),  A  is  an  Nt  x  Nj  matrix  whose  ij th  component  specifies 
whether  Ej.  has  reached  the  coercive  value  ECi.  To  specify  A,  we  define  the  matrices 


'  -l  •• 

•  -l 

i 

•  •  i  ' 

1 _ 

1 - 

^ init  — 
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Ni  x  Nj 

ii 

.  Ec*i 

EcNi  _ 

Ek  +  Ej  i  •  •  •  Ek  +  EIn  _ 
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Ek  +  Ej1 


Ek  +  Ejn . 
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Ni  x  Nj 


and  weight  vectors 


viui (EC1),  •••  ,vNivi(ECN.) 


lxNi 


WT  = 


wiu2{Eh), 


,wNjv2(EiN.) 


1  xJVj 


where  Ek  =  E(tk)  is  the  fcth  value  of  the  input  field.  The  points  ECi  and  Ej  in  the  definitions  are 
determined  by  the  quadrature  rule  being  employed  on  intervals  [0,  ECmax]  and  [Eimin,  Ejmax\  chosen 
according  to  the  physical  decay  conditions  (12)  —  i.e. ,  the  densities  v\  and  v2  are  negligible  outside 
these  regions. 

The  polarization  Pk  ~  P(Ek)  is  specified  by  Algorithm  1.  In  this  algorithm,  .*  indicates  compo¬ 
nentwise  matrix  multiplication  and  sgn  denotes  the  signurn  function.  The  first  step  in  the  for-loop 
updates  A  by  incorporating  the  status  of  the  previous  coercive  field  switch. 

Depending  on  the  methods  used  for  programming,  the  use  of  Algorithm  1  rather  than  utilizing 
conditional  if-then  constructs  can  reduce  runtimes  by  factors  in  excess  of  100  so  that  full  multiloop 
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model  simulations  run  in  the  order  of  seconds  on  a  workstation.  This  level  of  efficiency  is  necessary  to 
achieve  real-time  implementation  of  control  algorithms  utilizing  the  model.  Finally,  we  note  that  the 
algorithms  reported  in  [2]  for  the  kernel  (10)  which  incorporates  thermal  activation  are  on  the  order  of 
2-3  times  slower  than  Algorithm  1.  While  this  produces  analogous  reductions  in  the  speed  of  inverse 
algorithms,  the  resulting  model  is  still  sufficiently  efficient  to  facilitate  real-time  implementation. 

Algorithm  1. 

A prev  —  ^init 
for  k  =  1  :  Nk 

A  =  sgn(£fc  T  £c.  *  Aprev) 

P  =  v/k  +  PrA 
Pk  =  VTPW 

Aprev  —  A 

end 

3  Inverse  Relation  Between  Polarization  and  Field 

The  model  (11)  and  discretized  model  (13)  quantify  the  relation  between  input  fields  and  the  polar¬ 
ization  generated  in  ferroelectric  materials.  To  construct  an  inverse  filter  of  the  type  illustrated  in 
Figure  3,  it  is  necessary  to  quantify  the  inverse  P-E  relation.  We  summarize  here  an  algorithm  which 
accomplishes  this  when  the  E-P  relation  is  assumed  monotone  and  illustrate  the  filtering  process 
depicted  in  Figure  3  through  a  numerical  example.  Extension  of  the  algorithm  to  accommodate  the 
non-nronotone  field-displacement  behavior  shown  in  Figure  2  is  addressed  in  Section  5. 

The  first  step  in  the  construction  of  an  inverse  filter  involves  the  determination  of  an  initial 
(Eq.  P(j)  value.  This  is  typically  done  with  Eq  =  0  so  that  Pq  =  APr  is  the  positive  or  negative 
renranence  value  or  Pq  =  0  for  depoled  materials.  The  values  of  APr  can  be  computed  using  the 
forward  model  (11)  or  (13)  by  taking  one  step  to  AEmax  chosen  sufficiently  large  that  all  dipoles 
have  switched  and  then  stepping  back  to  E  =  0.  For  a  specified  value  of  P,  nronotonicity  in 
the  E-P  relation  is  exploited  and  the  forward  model  is  subsequently  advanced  until  the  prescribed 
polarization  is  crossed.  Interpolation  is  then  used  to  specify  a  final  field  value  corresponding  to  the 
prescribed  polarization.  This  process  is  outlined  in  Algorithm  2  where  specified  polarization  values 
are  designated  by  {Pk}  and  computed  values  by  {Pk}  for  k  =  1,  •  •  •  Nk. 

Algorithm  2. 

for  k  =  2  :  Nk 

Specify  Estep  >  0  as  fixed  or  adaptive 

rip  —  p*  _  p* 

al  ~  rk  rk- 1 

A E  =  dP  ■  Estep 
Etmp  —  Pk—  1  i  Pt.m.p  —  Pk—  1 
while  sgn (dP)  ■  (Pk  -  Ptmp)  >=  0 
Etmp  =  Etmp  t  A  E 

Ptmp  given  by  (13) —  e.g.,  as  implemented  in  Algorithm  1 

end 

Ek  given  by  linear  interpolation 


end 


The  flexibility  and  robustness  provided  by  the  inverse  Algorithm  2  are  illustrated  in  Figure  5. 
The  polarization  plotted  in  Figure  5(c)  is  employed  as  input  to  Algorithm  2  to  yield  the  P-E  relation 
plotted  in  Figure  5(a).  At  each  time  step,  the  resulting  field  value  is  then  employed  as  input  to  the 
forward  Algorithm  1  to  yield  the  E-P  curve  shown  in  Figure  5(b).  These  output  polarization  values 
Pout  are  compared  with  inputs  Pin  in  Figure  5(c)  and  the  absolute  errors  | Pin  —  Pout\  are  plotted  in 
Figure  5(d). 

From  these  results,  a  number  of  conclusions  can  be  drawn,  (i)  We  first  note  that  the  model 
and  its  inverse  provide  the  capability  for  characterizing  a  wide  range  of  symmetric  and  biased  minor 
loop  behavior  —  e.g.,  see  [21,23,30].  (ii)  The  composition  of  the  inverse  and  model  in  the  manner 
depicted  in  Figure  3  can  effectively  linearize  the  nonlinear  transducer  behavior  with  the  numerical 
accuracy  | pn  —  Pout\  limited  only  by  the  stepsize  dP.  Whereas  the  accuracy  in  a  physical  system  will 
be  diminished  due  to  modeling  error,  linearization  in  this  manner  can  significantly  improve  control 
authority  since  less  control  effort  is  focused  on  unmodeled  or  nonlinear  dynamics.  This  forms  the 
crux  of  various  linear  control  designs  [15, 16].  (iii)  Although  faster  implementation  algorithms  can  be 
constructed  for  the  inversion  process  [8] ,  the  algorithm  described  here  is  highly  robust  and  avoids  the 
potential  for  losing  track  of  the  memory  incorporated  in  the  model.  Furthermore,  the  use  of  adaptive 
stepsizes  A E  ensures  that  Algorithm  2  is  approximately  a  factor  of  two  slower  than  the  forward 
algorithm  which  is  reasonable  for  physical  implementation,  (iv)  Whereas  Algorithm  2  employs  the 
limiting  piecewise  linear  kernel  P  given  by  (4),  analogous  algorithms  have  been  developed  for  the 
more  general  kernel  (10)  which  incorporate  thermal  relaxation  and  additional  dynamic  effects  [8]. 


-2-1  0  1  2 
Electric  Field  (MV/m) 


(a) 


(b) 


Normalized  Time  (sec) 

(c) 


Normalized  Time  (sec) 

(d) 


Figure  5:  (a)  Inverse  relation  Pin-Eout  given  by  Algorithm  2.  (b)  Forward  relation  Eout-Pout  from 
Algorithm  1.  (c)  Comparison  between  P,;n  and  Pout •  (d)  Absolute  error  |Pm  —  Pout\  for  complete 
inversion  process  depicted  in  Figure  3. 
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4  Lumped  Model  for  the  Stacked  Actuator 


The  constitutive  relation  (15)  quantifies  the  electromechanical  behavior  of  piezoceramic  materials 
operating  below  the  coercive  stress  ac  where  ferroelastic  switching  commences.  In  this  section  we 
use  this  relation  to  construct  a  macroscopic  model  for  the  stacked  PZT  rod  employed  as  an  AFM 
stage  in  the  manner  depicted  in  Figure  1(a).  While  this  illustrates  macroscopic  model  development 
for  a  specific  application,  similar  principles  hold  for  other  transducer  designs  —  e.g.,  see  [21]. 

We  assume  that  the  stacked  actuator  or  rod  has  length  £,  cross-sectional  area  A,  density  p, 
Young’s  modulus  Y,  and  Kelvin- Voigt  damping  parameter  cd .  We  also  assume  that  the  end  at 
x  =  0  is  fixed  whereas  the  end  at  x  =  £  is  subjected  to  inertial,  elastic  and  damping  forces  associated 
with  the  stage  mechanisms.  Because  material  properties  and  forces  along  the  length  of  the  rod  are 
uniform,  we  consider  a  lumped  model  quantifying  the  displacement  u(t)  at  x  =  £.  The  validity 
of  the  lumped  ODE  model  as  compared  with  a  distributed  PDE  model  is  established  in  [24],  The 
geometry,  mass  mg,  stiffness  k ^  and  damping  mechanisms  eg  associated  with  the  end- forces  at  x  =  £ 
are  depicted  in  Figure  6. 

From  the  assumption  of  uniform  stresses  and  strains  through  the  length  of  the  rod,  it  follows  that 


e(t)  = 


u{t) 

£ 


in  the  stress  relation  (15).  Balancing  the  forces  a  A  for  the  rod  with  those  of  the  stage  mechanisms 
yields  the  lumped  model 

,d2u,  .  cnAdu ,  .  Y  A  ,  ,  d2u ,  .  du ,  .  ,  .  . 

M-^-(t)  +  u(t)  -  AaiP(E(t ))  -  Aa2P2{E(t))  =  -  ce~{t)  -  ku(t ) 

or,  equivalently, 

+  c  “77(f)  +  ku(t)  =  diP(E(t))  +  d2P2{E(t))  (16) 

dtz  at 

where 

cd  A  Y  A  „  . 

m  =  pA  +  mu  ,  c  =  — - - h  eg  ,  k  =  — — \-  kc  ,  a\  =  Aa\  ,  a2  =  Aa2 

and  the  initial  conditions  are  u(0)  =  uq  and  ^(0)  =  u  1 .  The  polarization  P  is  specified  by  the  model 
(11)  or  discretized  model  (13). 

The  model  can  also  be  written  as  the  first-order  system 


u(t)  =  Au(t )  +V(E(t)) 
u(  0)  =  u0 


(17) 


u(t) 

_ 


CE- 


x  =  0 


x=Z 


Figure  6:  Rod  geometry  used  when  modeling  the  stacked  PZT  actuator  employed  in  the  AFM  stage 
depicted  in  Figure  1(a). 
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where  u(t )  =  [u(t) ,  u(t)]T ,  u{ 0)  =  [uq,ui]t  and 

,  V(E(t))  =  ^[a1P(E(t))+d2P2(E(t))\  J  . 

This  formulation  proves  advantageous  in  the  next  section  when  establishing  notation  used  in  the 
construction  of  the  inverse  algorithm. 

The  accuracy  of  the  framework  is  illustrated  in  Figure  7  where  the  lumped  model  (16)  with 
P  specified  by  (13)  is  used  to  characterize  the  frequency-dependent  dynamics  of  the  PZT  stacked 
actuator  employed  in  the  AFM  stage  depicted  in  Figure  1(a).  When  constructing  the  polarization 
model,  the  general  densities  v\  and  were  identified  using  the  least  squares  techniques  detailed 
in  [23].  It  is  noted  that  the  combined  model  quantifies  both  the  hysteresis  and  dynamic  effects 
observed  as  frequencies  are  increased.  Additional  details  regarding  the  construction  and  validation 
of  the  stacked  actuator  model  for  the  AFM  stage  are  provided  in  [8,24]  whereas  additional  examples 
demonstrating  properties  of  the  model  for  characterizing  hysteresis  in  various  PZT  compounds  can 
be  found  in  [21-23,30]. 


A  = 


0  1 

—k/m  —c/m 


(a)  (b)  (c) 

Figure  7:  Characterization  of  AFM  field-displacement  behavior  with  sample  rates  of  (a)  0.279  Hz, 
(b)  5.58  Hz  and  (c)  27.9  Hz. 


5  Inverse  Relation  Between  Displacements  and  Fields 


The  inversion  algorithm  summarized  in  Section  3  relies  on  the  monotonicity  of  the  E-P  relation.  As 
illustrated  in  Figure  2,  this  property  is  not  shared  by  the  E-u  relation  as  frequencies  are  increased 
so  we  develop  here  an  extended  inversion  algorithm  which  incorporates  this  non-nronotone  behavior. 
The  crux  of  the  modification  focuses  on  the  accommodation  of  dynamic  effects  in  the  E-u  behavior. 

To  establish  notation  used  when  quantifying  dynamic  effects,  we  employ  modified  semigroup 
notation  to  define  solution  values 


Uk+i  =  u(tk+1,tk,E,  uk)  =  CeA(-tk+1  tk)uk  +  C  eA(tfc+1  s^V(E(s))ds 

Jtk 

i 

Uk+1  =  u(tk+i,tk,Ek,uk)  =  CeA^tk+1~tk)uk  +  c  eA{tk+1~s)V{Ek)ds 


'tk 


(18) 


where  C  =  [1,  0]  and  Ek  =  E{tk).  Hence  uk+\  is  the  solution  to  (16)  or  (17)  with  the  electromagnetic 
force  applied  throughout  the  time  interval  [tk,tk+ i]  whereas  uk+\  denotes  the  displacement  of  the 
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rod  at  x  =  £  with  the  force  maintained  constant  at  the  kth  value  throughout  the  interval  [tk,tk+ 1]  as 
depicted  in  Figure  8.  The  definition  of  Uk  is  similar  to  that  of  uk+ 1- 
The  upper  and  lower  hysteresis  curves  are  denoted  by 

r ,  ^  ..dE  1 

C1  =  {(-B-“)I  It  -°j 

W  =  <o}^ 

Due  to  dynamic  (inertial)  effects,  it  is  observed  that 


Uk+l  —  OH  C-L 

Uk+ 1  >  Uk+i  on  cjj 


(19) 


as  depicted  in  Figure  8(b)  and  (c).  To  motivate  the  relation  on  cl,  we  note  that  (18)  yields 

|_i  p  i 

Uk+ 1  -  uk+ 1  =  C  /  eA^+i-s)  V(E(s))  -  V(Ek)  ds.  (20) 

Jtk  1 

We  now  establish  that  the  right  hand  side  of  (20)  is  nonnegative  given  that  the  monotonicity 
of  the  E-P  relation  implies  that  V(E(s))  —  T{Ek)  >  0  for  s  >  tk.  We  first  note  that  the  Cayley- 
Hamilton  theorem  dictates  that  eA^tk+1~s^  =  olqI  +  ol\A.  Furthermore,  it  follows  from  the  definitions 
of  C  and  V(E(t))  that  only  the  (1,2)  entry  of  eA<'tk+'~si  contributes  to  the  right  hand  side  of  (20). 
Moreover,  it  follows  from  the  definition  of  A  that  the  (1,2)  entry  of  ^Mtk+i-s)  jg  sjmpiy  ai. 

To  determine  a\ ,  we  note  that  there  are  two  possibilities  for  the  eigenvalues  of  A:  (i)  both  are 
real,  distinct,  and  negative,  and  (ii)  they  are  a  conjugate  pair  with  negative  real  part.  We  consider 
the  first  case  where  the  eigenvalues  satisfy  A2  <  Ai  <  0.  It  follows  from  the  Cayley-Hamilton  theorem 
that 

e\i(tk+i-s)  _  aQ  _|_  ai\1 
eA2(ife+i-s)  _  aQ  _(_  ai\2 


(b)  (c) 

Figure  8:  Non-monotone  behavior  of  the  field- displacement  relation  measured  at  27.9  Hz  as  shown 
in  Figure  2(c).  Solution  values  uk+ 1  and  uk+ 1  respectively  due  to  field  inputs  Ek  and  E(t)  on  the 
(b)  lower  loop  ck.  and  (c)  upper  loop  cjj  of  the  curve. 
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so  that 


Ot\ 


g^l(tk+l  _  g^2(^fc+l 

Ai  —  A2 


>  0. 


Similar  analysis  holds  for  the  second  case.  Thus  the  integrand  CeA^tk+1~s^[V(E(s))  —  V(Ek)\  of  (20) 
is  nonnegative  and  consequently  so  is  the  integral.  As  a  result,  uk+\  —  uk+\  >  0  or  equivalently 
uk+ 1  <  Uk- 1-1.  The  argument  for  cjj  is  analogous. 

In  the  inverse  algorithm  used  to  specify  the  u-E  relation  given  data  values  {u£},  k  =  1,  •  •  •  ,Nk, 
the  inequalities  (19)  are  applied  to  either  the  exact  or  discretized  solution  of  (16)  to  determine  the 
appropriate  sign  of  A E  when  implementing  the  polarization  component  (13)  of  the  model.  The 
notation  utmp  and  utmp  designate  either  the  exact  or  approximate  solutions  to  (13)  having  the 
interpretation  specified  in  (18).  The  resulting  inversion  process  is  outlined  in  Algorithm  3.  Attributes 
of  the  algorithm  are  illustrated  in  the  next  section  in  the  context  of  experimental  open  loop  control 
implement  ation . 


Algorithm  3. 

for  k  =  2  :  Nk 

Specify  Estep  >  0  as  fixed  or  adaptive 
Specify  At 

A E  =  sgn (u*k+1  -  uk+ 1)  •  Estep 
Etmp  —  Ek_  1  ,  Pt.mp  —  Pk—1  1  ttmp  —  T—  1 
while  sgn(u*k+1  -  uk+ 1)  •  A E  >=  0 
ttmp  =  ttmp  T  At 

compute  utmp  as  true  or  approximate  solution  to  (16) 
Etmp  =  Etmp  T  A  E 

Ptmp  given  by  (13) —  e.g.,  as  implemented  in  Algorithm  1 

end 

Ek  given  by  linear  interpolation 


6  Open  Loop  Control  Implementation 

To  illustrate  the  effect  of  filters  employing  the  inverse  model  developed  in  Section  5  on  open  loop 
tracking  performance,  we  summarize  experiments  conducted  at  0.279  Hz  and  27.9  Hz.  The  trajectory 
to  be  tracked  consisted  of  triangle  waves  having  amplitudes  of  40.56  //m  and  27.04  /xnr  as  shown  in 
Figures  9-11. 

To  specify  the  model,  parameters  in  the  polarization  model  (13)  and  lumped  rod  model  (16) 
were  estimated  through  a  least  squares  fit  to  field-displacement  data  collected  at  0.279  Hz,  5.58  Hz 
and  27.9  Hz  using  the  techniques  detailed  in  [23].  This  yielded  model  fits  similar  to  those  shown 
in  Figure  7.  The  model  with  these  parameters  was  then  used  to  construct  an  inverse  filter  using 
Algorithm  3  of  Section  5.  In  a  series  of  experiments,  the  specified  trajectories  were  input  to  the  filter 
and  the  resulting  field  was  applied  to  the  AFM  stage.  To  provide  a  metric  for  comparison,  a  second 
input  field  for  each  case  was  determined  through  a  linear  scaling  of  the  field-displacement  relation. 
This  linear  filter  accommodated  the  scaling  difference  between  inputs  and  outputs  but  neglected 
inherent  hysteresis  and  constitutive  nonlinearities. 
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Figure  9:  Tracking  performance  utilizing  the  model-based  inverse  filter  of  Section  5  and  a  linear 
filter  at  0.279  Hz  and  amplitude  of  40.56  gm.  (a)  Specified  trajectory  and  tracking  provided  by  the 
inverse  and  linear  filters,  and  (b)  errors  obtained  with  the  two  filters),  (c)  Field- displacement  data 
measured  before  and  after  the  open  loop  control  experiment. 


(a)  (b) 

Figure  10:  Tracking  performance  utilizing  the  model-based  inverse  filter  of  Section  5  and  a  linear 
filter  at  0.279  Hz  and  amplitude  of  27.04  gm.  (a)  Specified  trajectory  and  tracking  provided  by  the 
inverse  and  linear  filters,  and  (b)  errors  obtained  with  the  two  filters). 
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Figure  11:  Tracking  performance  utilizing  the  model-based  inverse  filter  of  Section  5  and  a  linear 
filter  at  27.9  Hz  and  amplitude  of  27.04  ^m.  (a)  Specified  trajectory  and  tracking  provided  by  the 
inverse  and  linear  filters,  and  (b)  errors  obtained  with  the  two  filters),  (c)  Field- displacement  data 
measured  before  and  after  the  open  loop  control  experiment. 


The  specified  and  achieved  trajectories,  errors,  and  selected  hysteresis  plots  for  the  combinations 
(i)  0.279  Hz,  40.56  /jm  amplitude,  (ii)  0.279  Hz,  27.04  /jrn  amplitude,  and  (iii)  27.9  Hz,  27.04  //m 
amplitude  are  plotted  in  Figures  9-11.  At  0.279  Hz,  the  filtered  design  provides  only  marginally 
improved  accuracy  due  to  the  low  degree  of  hysteresis.  For  the  more  hysteretic  response  at  27.9  Hz, 
however,  the  inverse  filter  provides  a  significant  increase  in  accuracy  and  yields  errors  that  are 
approximately  a  factor  of  10  smaller  than  the  linearly  scaled  case.  This  illustrates  the  advantage 
of  incorporating  the  frequency-dependent  model  inverse  in  the  control  design.  We  note  that  the 
same  model  parameters  were  employed  in  each  case  thus  illustrating  the  capability  of  the  model  and 
model-based  inverse  to  compensate  for  the  frequency-dependent  hysteresis. 

The  primary  source  of  errors  in  the  filtered  design  is  variability  between  experiments  as  illustrated 
by  the  variation  in  the  hysteresis  plots  measured  at  the  two  frequencies  before  and  after  the  open 
loop  control  experiments.  This  is  hypothesized  to  be  due  to  variations  in  the  true  applied  voltage 
and  illustrates  one  reason  feedback  is  necessary  in  final  control  designs. 

7  Concluding  Remarks 

This  paper  addresses  the  development,  implementation  and  experimental  validation  of  a  model-based 
inverse  filter  to  accommodate  hysteresis  and  constitutive  nonlinear  inherent  to  the  field-polarization 
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and  field-displacement  behavior  of  ferroelectric  materials.  The  nonlinear  hysteresis  effects  are  quan¬ 
tified  using  a  previously  developed  framework  consisting  of  energy  relations  at  the  lattice  level  in 
combination  with  stochastic  homogenization  techniques  to  provide  low-order  macroscopic  constitu¬ 
tive  relations.  The  development  of  a  lumped  transducer  model  based  on  these  constitutive  relations 
is  illustrated  in  the  context  of  a  PZT  stage  for  an  AFM.  The  inverse  displacement-field  model  exploits 
monotonicity  in  the  E-P  relation,  the  efficiency  of  forward  E-P  algorithms,  and  dynamic  properties 
of  the  transducer  model. 

To  illustrate  attributes  of  the  inverse  displacement-field  algorithm,  it  was  employed  as  a  filter 
in  open  loop  tracking  experiments  with  an  AFM  stage.  These  experiments  illustrate  that  at  low 
frequencies,  where  hysteresis  and  constitutive  nonlinearities  are  minimal,  incorporation  of  the  in¬ 
verse  filter  provides  only  marginal  improvement  in  tracking  accuracy  as  compared  with  a  linear 
filter.  However,  at  higher  frequencies  where  hysteresis  becomes  significant,  the  inverse  filter  yields 
an  approximately  tenfold  improvement  in  accuracy  compared  with  the  linear  filter  thus  maintaining 
tracking  accuracy  even  though  the  transducer  is  operating  in  highly  hysteretic  and  nonlinear  regimes. 
Linearization  of  the  electromechanical  behavior  in  this  manner  reduces  the  degree  to  which  feedback 
mechanisms  must  expend  energy  linearizing  the  transducer  response  and  increases  control  authority 
for  stabilization  or  tracking. 

We  note  that  a  significant  advantage  of  the  energy-based  model  is  the  fact  that  it  provides  a  unified 
framework  for  characterizing  hysteresis  and  constitutive  nonlinearities  in  ferroelectric,  ferromagnetic 
and  ferroelastic  (e.g.,  SMA)  compounds  [28,29].  One  facet  of  present  investigations  focuses  on 
the  extension  and  implementation  of  the  inverse  filtering  techniques  for  the  latter  two  classes  of 
compounds. 

Present  investigations  are  also  focused  on  the  development  of  robust  feedback  control  designs 
which  employ  the  inverse  filters  to  linearize  transducer  dynamics.  Initial  investigations  focused  on 
the  numerical  implementation  of  these  energy-based  inverse  filters  for  magnetic  transducers  will  be 
reported  in  [17]  and  the  experimental  implementation  of  feedback  designs  exploiting  the  filters  are 
under  present  investigation. 
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